NASA Contractor Report 185676 


//. ■■/- 
P t /(. 

Reference Equations of Motion for 
Automatic Rendezvous and Capture 


David M. Henderson 

TRW Systems Services Company, Inc. 

Houston, Texas 

(TRW Report Number 91:J431.1-182) 


Prepared for 

Lyndon B. Johnson Space Flight Center 
under contract NAS9-17900 


1992 


EQUAr7«; , 22 6 :‘» IO «FMF'*e 

***««<*« ANo'cIpfuaT 


AUTOMATK 
C TRW ) 


N93-29652 


Unci as 


G3 /13 0176443 




NASA Contractor Report 185676 


Reference Equations of Motion for 
Automatic Rendezvous and Capture 


David M. Henderson 


Contract NAS9-17900 


April 1992 




TABLE OF CONTENTS 


1.0 INTRODUCTION 1 

2 . 0 TRANSLATION DYNAMICS - MOTION OF THE CENTER 

OF MASS 3 

3 . 0 ROTATIONAL DYNAMICS - VEHICLE ANGULAR MOTION 

ABOUT THE CENTER OF MASS 9 

4.0 SIX DEGREES OF FREEDOM EQUATIONS OF MOTION 11 

5.0 TRANSLATION CONTROL PARAMETERS 14 

6.0 ROTATION CONTROL PARAMETERS 21 

7.0 CONCLUSIONS AND RECOMMENDATIONS 2 5 

8.0 REFERENCES 26 

APPENDIX A SUGGESTED INTEGRABLE PARAMETER LIST 2 7 

APPENDIX B GRAVITATIONAL ACCELERATION 2 9 

APPENDIX C PARTIAL DERIVATIVES OF THE GRAVITATIONAL ACCELERATION VECTOR 3 2 


Mi 


\ 

m i s , , w 




REFERENCE EQUATIONS OF MOTION FOR 
AUTOMATIC RENDEZVOUS AND CAPTURE 

David M. Henderson 

TRW/Houston 


1.0 INTRODUCTION 

The analysis presented in this paper defines the reference 
coordinate frames and control parameters necessary to model the 
relative motion and attitude of spacecraft in close proximity with 
another space system during the Automatic Rendezvous and Capture 
phase of an on - orbit operation. Translation motion equations 
relative to the geodesic LVLH frame are developed yielding the 
Clohessy-Wilt shire differential equations of motion. 

Perturbations of the gravity partial derivative matrix due to a 
4x4 gravity model show the expected deviations from the Clohessy - 
Wiltshire approximation. For the docking case, however, as the 
relative separation distances become small, both relative 
gravitational and inertial accelerations can be omitted from the 
translation equations of motion. 

Euler’s equation for the rotational rate accelerations of a rigid 
body are developed with expected torque effects indicated. The 
actual torque computations are problematic and these specific 
models are not developed as part of this work. 

The relative docking port target position vector and the attitude 
control matrix are defined based upon an arbitrary spacecraft 
design. These translation and rotation control parameters could 
be used to drive the error signal input to the vehicle flight 
control system. Measurements for these control parameters would 
become the bases for an autopilot or FCS design for a specific 
spacecraft . 

Cartesian analysis is used throughout the paper for the clarity 
and compactness in representing component equations . The use of 
Cartesian tensor notation makes the task of repeated 
differentiation of transformation equations considerably easier 
and more understandable. For instance, the position vector in 
inertial space is given by 


x = a x + X 


( 1 . 1 ) 


where ax is a transformation matrix times the vector x. Equation 
1.1 defines both rotation and translation in the inertial frame. 
Differentiating Equation 1.1 for velocity is simply, 
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v = au + V 


( 1 . 2 ) 


where 


v = W X + v 


(1.3) 


and W is the skew-symmetry angular velocity matrix. The vector v 
is recognized from vector analysis [see Reference 4.0, Page 208] 
as the operator. 


aiX 


+ 


d > 

dt j 


(1.4) 


However, in Equation 1.2 the inertial coordinates are preserved 
and no ambiguities arise as to which coordinate system the vectors 
components are represented. 
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2.0 TRANSLATIONAL DYNAMIC MOTION OF CENTER OF MASS 


Let (o) be the origin of an imaginary coordinate frame attached to 
a mass in orbital motion about the planet^s center . farther 
assume that there is no contact force acting upon this mass s 
lilt i t moves along the geodesic arc in the inertial frame whose 
origin coincides with the center of mass of the planet. Also 
attached to the origin (o) is the Local-Vertical-Local-Horizontai 
(LVLH) coordinate reference frame. For circular oroital motion 
the LVLH system rotates downward about its y-axis at the orbital 

rate . 



The Cartesian equation, 


— a^x. 


+ X r 


( 2 . 1 ) 


Represents the 
center located 
transformation 


components of the inertial vector to any coordinate 
at (n) as seen from (o) , where a Q is the 
matrix LVLH to inertial and 
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3-o^on X 0 


(2.2) 


Xn = 


represents vectors to the center of mass m„ at (n) or 


a Q 


k on 


X 


o • 


(2.3) 


This is their relative position relationship at any time t. 
Differentiating Equation 2.3 for velocity 


a G ^on — V n — V 0 


(2.4) 


where 


^on = W X on + V on 


and 


0 

“®3 

©2 


0 

-©! 

— fi>2 

©1 

0 


is the skew-symmetric angular velocity matrix. Solving for 
relative velocity in the LVLH frame. 


?o„ = (a G ) T (V n -V c ) - WX on . 
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Differentiating Equation 2.4 for acceleration 


^o^on V 0 


(2.5) 


where 


a on = (WW + W) X on + 2W ? on + V 0 


and the relative LVLH acceleration vector as seen from (o) is 


^on = (a 0 ) T (V n -V 0 ) - (WW + W) X on - 2 W V on . 


( 2 . 6 ) 


The acceleration vectors are 


V n = G(X n ) + A n 


and 


V 0 = G(Xo) + Ao 


(2.7) 


But Aq = 0, since we have defined (o) to be an inertial frame 
moving along the geodesic arc about the planet center . Therefore* 


V n - V G = G(X n ) - G(X 0 ) + A n 


( 2 . 8 ) 


and Equation 2.5 becomes 


a G a on = G(X n ) - G(X 0 ) + A n . 


(2.9) 


If X n and X 0 are very close to one another. The difference in 
gravitation accelerations is simply a linear variation in 
acceleration about the origin (o) , 
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SG(X 0 ) = G(X n ) - G(X 0 ) 


and 


6G SG 

5G(X 0 ) = SX 0 = (X n -X 0 ) = P x (X n -X 0 ) , (2.10) 

SX 0 6X 0 


therefore. Equation 2.9 becomes 


a oSon Px (X n -X Q ) + A n . 


( 2 . 11 ) 


From Equation 2.3 we have 


a oSon Px ( a oXon) ^ A n 


( 2 . 12 ) 


and solving for the acceleration in the LVLH frame. 


^on — [ (a Q ) T Px a e ] ^on + (a G ) T A n 


(2.13) 


a 


on 


Px X on + A n 


(2.14) 


where 


A n — ( a o) T A n and Px — t( a o) T Px a ol • 
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The P x matrix is defined from the similarity transformation for 
matrices as described in Reference 5.0, Pages, 317 and 318, then 
Equation 2.14 becomes. 


Von = (Px-WW-W) X on - 2W V on + A n 


(2.15) 


Equation 2.15 is the equation for motion for mass m in the LVLH 
coordinates system attached to the geodesic point (o) in Figure 
1.0. The P x matrices are used here for simplicity and also 
contribute to understanding in the analysis, however, it should be 
pointed out that the ^ x matrices are non-linear and their region of 
numeric stability about the geodesic origin is small. Hence, for 
very accurate relative motion computations the difference 
equations of Equation 2.9 become a more reliable relationship on 
which to base computational algorithms. Equation 2.15 can be 
represented by the sum of two vectors. 


Von = A ln + A n, 


(2.16) 


where Ai n are the combined gravitational and inertial accelerations 
effecting the relative motion of the center of mass . 

The P x matrix reduces to (see Appendix C) 



'l 

0 

0 ^ 


r l 

0 

0> 

P x = -HR' 3 

0 

1 

0 

V 

1 

II 

0 

1 

0 



0 

"2 j 


,0 

0 

"2; 


In a spherical gravitational field, i.e., where all gravity 
harmonics are assumed zero, 
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r o 

0 

-®o 

\ 


'o 

0 

> 

-Wo 


r -wg 

0 

0 ' 

ww = 

0 

0 

0 



0 

0 

0 

= 

0 

0 

0 


v<°o 

0 

0 

) 


^<>o 

0 

0 , 


.0 

0 

-wg > 


(2.18) 


and for uniform circular 


orbit motion 


W 


0 . 


(2.19) 


However, in the "real world" small oscillatory motions do occur 
about the LVLH axes due to the gravity harmonics. Then Equation 
2.15 reduces to 


y x *= 2®o V z + A x = A Ix + A x 

Vy = -<»g y + Ay = A Iy + Ay 

V z = 3cog 2 - 2o 0 V x + A z = A Iz + A z 


( 2 . 20 ) 


which are recognized as the Clohessy-Wiltshire approximation for 
the differential equations of relative motion [as described in 
Reference 1.0]. The Aj in Equation 2.20 are the combined relative 
gravitational and inertial acceleration effects in the LVLH 
system. 
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3 . 0 ROTATIONAL DYNAMICS 


The vehicle angular momentum in body axis coordinates transforms 
into the inertial frame via 


and L for any rigid body is given by Ia> where I is the body axis 
moment of inertia matrix and ro is the body axis rotation rate 
vector. Equation 3.1 becomes 


Using Newton's Law; 


dL 

— = N. 
dt 


(3.3) 


Equation (3.2) becomes 


UlJ _ . _ . _ . . 

— = aWI«> + alffi = N = aN (3.4) 

dt 


where N and N are the applied torques in the inertial and body 
axis coordinate frames, respectively. Equation 3.4 yields Euler's 
Equation of motion for a rigid body, 

n = wra + Is, (3.5) 


and solving for 5, we have, 
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© = J (N-WI©) 


(3.6) 


where J = I 1 , is the inverse of the body axis moment of inertia 
matrix. The applied torques are composed of 


N = N G + N T + N A + N P 


(3.7) 


gravity, thrusting, aerodynamic and plume impingement torques, 
respectively, and Equation 3.6 becomes 


© 


J 


[N g + N x + N a + Np + Nj] 


(3.8) 


where the inertial torques are Nj = -WI5, which results due to 
rotational motion of the vehicle itself. Integrals of this 
equation [see Reference 3.0] provide the body axis rotation rate 
vector as a function of time. 


10 


4.0 SIX DEGREES-OF-FREEDOM EQUATIONS OF MOTION 
From Equation (2.20), 

9o 

v n = — X F n + A In 
W n 

for translational dynamics and Euler's equation 
(b = Jn (N n W n I n «>n ) 
for rotation dynamics, where; 


(4.1) 


(4.2) 


— X F n " Summation of contact forces acting on vehicle n 
W n (LVLH Ref) 


Ain 


In 

Jn 


W n 


<b 


“n 


Inertial and gravitational accelerations of CM 
of vehicle n (LVLH Ref) 

~ Moment of inertia matrix of vehicle n 

~ Inverse of moment of inertia matrix of vehicle n 

“ Body axis torques about the CM of vehicle n 

Skew symmetric angular velocity matrix 
(vehicle n) 

“ Body axis rotation rate vector (vehicle n) 

" Body axis rotation rate acceleration vector for 
vehicle n. 



For example an integrable set of differential equations could be 
represented as follows. 


• 9o 

Vi = (~Cdi q Si Qreli+T!) + An 

Wi 


(4.3) 


• 9o 

V 2 = (~Cd2 <5 ^2 QreL 2 ) + A12 . 

W 2 


(4.4) 


Assume that vehicle one is active (with thrust) and that vehicle 
two is passive (no thrusting) and 


• 9o 

^1 = “®i9o 3 Ureli + Ti + An 

Wi 


V2 ~ -B 2 g 0 Q UREL2 + Aj 2 , 


(4.5) 


where u REL is a unit vector in the LVLH geodesic system and is in 
the direction of the motion of the geodesic frame and, 


Bi is the ballistic coefficient 


r C Di Si ^ 
^ Wi ; 


q is the fluid dynamic pressure, q = (1/2 p (V R ) 2 ) where, 


V R " relative velocity magnitude 


p ” fluid density 
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and the u REL vector in the LVLH reference frame is approximately. 


1.0 
0.0 
0.0 

With the above example in mind. Appendix A shows a suggested 
integrable parameter list for a time - wise step solution to the 
problem. More precise equations of motion can be defined by 
computing the driving values for F n and N n in Equations 4.1 and 
4.2, respectively, however, the same 6-DOF solution procedures 
will apply. 
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5.0 TRANSLATION CONTROL PARAMETERS 

As in Figure 2-1, the geodesic center is the motion of an 
imaginary point falling along the local gravity geodesic arc. No 
drag or other contact forces are acting at the geodesic center. 
The rotating LVLH coordinate axes are attached to the geodesic 
center. The position vector in the inertial coordinate system to 
a point in space can be computed from each of the coordinate 
systems shown in the following figure. 



X 

Figure 5-1 
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The Cartesian equations relating the coordinates as shown in 
Figure 5.1 can be written as follows; 


x 0 = a 0 x 0 + X 0 (5.1) 

from the geodesic frame and likewise, 

xi = ai*! + Xi (5.2) 

is from the center of mass of vehicle one. From coordinates xi on 
vehicle one, 

Xi = c a X x + x lf (5.3) 

is from the center of the docking port located on vehicle one. 
Similarly for vehicle two, 

X 2 = 32^2 + ^2 (5.4) 

from the center of mass of vehicle two 

X 2 = C 2 X 2 + X 2 (5.5) 

and from the center of the docking port located on vehicle two. 
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The following definitions apply to Equations 5.1 through 5.5 


a Q “ The LVLH to inertial transformation matrix attached to 
the geodesic coordinate center {a time varying matrix) 


ai " i th body axis to inertial transformation matrix 
(a varying matrix) 


Xi ” The inertial coordinates of the i th coordinate center 
(a dynamic vector in motion about the planet center) 

Ci " i th docking port to body axis transform matrix 
(a constant matrix) 

Xi ~ Body axis coordinates of the i th body docking port 
coordinate center (a constant vector) . 


The target docking port inertial position is 


Xt — &2 ^2 ^ X 2 


x T = a 2 (c 2 X 2 +X 2 ) + X 2 


(5.6) 


but 


x 2 = 0, 


since the target is the docking port coordinate center and 
therefore Equation 5.6 reduces to 


xt = a 2 X 2 + X 2 , 


(5.7) 


the target location in inertial coordinate. Hence, the position 
relative geometry in the inertial coordinate reference frame 
becomes 


1 6 


a 2 X 2 + X 2 = a.i (Ci*i+Xi) + Xi (5.8) 


from Equation 5.1 


X! = a 0 X 0l + X 0 and X 2 = a 0 x o2 + X Q (5.9) 


and substituting these into Equation 5.8 


a 2 X 2 + a 0 x o2 = a x (Ci Xj+Xi) + a 0 x 0l (5.10) 


where the X o2 and X o2 are the position vectors to the vehicle one 
and two center of mass measured in the LVLH geodesic coordinate 
frame of reference . 

By multiplying Equation 5.10 by the transformation matrix (a 0 ) T , 
the equation is transformed to the LVLH frame, 

(a c ) T 32 X 2 + X 0 2 = (a c ) T ai (ciX^+Xi) + x 0 2 « (5.11) 

Now let bi = (a 0 )Tai and b 2 = (a 0 ) T a 2 , be the transformation matrices 
for the body axes to the LVLH reference frame. Then Equation 5.11 
becomes 

(X 0 2 — X 0 i) = bi (CiXi+Xi) - b2X 2/ (5.12) 

the relative position vector between the two vehicles in the LVLH 
frame and solving for 

Xi = (Ci)T[ (bi)T{ (X 02 -X 0l )+b2X 2 }- XJ . (5.13) 
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We have the position control vector as seen from the vehicle one 
clocking port coordinate center. The GN&C systems should drive 
this vector to zero for docking contact. Written in body axis 
coordinates we have. 


Xl = Cl*! = (bx)T (b 2 X 2 + (x o 2 -x o1 ) ) - X lf 


(5.14) 


the control vector as seen from the active vehicle (vehicle one) . 
Hence, a candidate commanded velocity in the vehicle frame 
becomes. 


v c = K r X x 


(5.15) 


where the gain. 


Kr 



V4t ; 


can be applied in an Flight Control System (FCS) with appropriate 
limiters for proper closing rates . For the velocity of the target 
docking port we differentiate Equation 5.9 


a 2 W 2 X 2 + a 0 W 0 x 02 + a 0 v o2 = a^ (Ci^x+Xx) 


+ a^i?! + a 0 W 0 x 0l + a 0 v 0l . 


(5.16) 


Again by multiplying by (a 0 ) T the equation is reduced to LVLH 
coordinates. 
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(a Q ) T a 2 W 2 X 2 + W 0 x o2 + v o2 = (a c ) T aiWi (Cx^i+Xi) 


+ (a^TaxCi^i + W 0 x 0l + v o1 


(5.17) 


and the relative velocity difference in the LVLH frame is. 


(^ 02 -^ 01 ) = bxWx (cxXx+Xx) + biCi^x - b 2 W 2 X 2 - W 0 (x o 2 -x 0 i) (5.18) 


and the docking port relative velocity of the docking port 
coordinate center on vehicle two (the target) is 


= (cx)T[ (bx)T{ (v 02 -y 0 i)+w 0 (x 0 2 -x 0l )+b 2 w 2 x 2 }-w 1 (c 1 x 1 +x 1 ) ] . (5.19) 


The control velocity vector in vehicle body axis coordinates 
becomes, 


?x = ci^x = (bi) T [ (v 02 -y 0l ) +W Q (x o2 -x o1 ) +b 2 w 2 x 2 ] -Wx (CxSx+Xx) . (5.20) 


Using Equations 5.15 and 5.20, the AV command for translation 
control becomes. 


av c = tr c - Vx, 


(5.21) 


and when applied with appropriate limiters or closing gates, the 
algorithm can be used to the control closing conditions for safe 
docking. The acceleration of the control vector is found by 
differentiation of Equation 5.16 and becomes 



= (Ci)T [ (bi)T b 2 (W 2 W 2 +W 2 )X 2 +(b 1 )T{ (W 0 W 0 +W 0 ) (X o 2 -Xol> 


+ 2 


(Wo) (v o 2 -v ol ) + (v o 2 -v ol ) 


} 


- (WiWi+Wx) (CiXi+Xi)-2Wc!^i] 


(5.22) 


and in body axes coordinates 


(ci^i) - (bi)T b 2 (W 2 W 2 +W 2 ) X 2 + (bi)T { (W 0 W 0 +W 0 ) (x o 2 -x 0 i) 


+ 2 


W 0 (^o2-^ol) + (^o2-^ol) > 


-(WiWi+Wi) (dXi+Xi)^ WiCi^x. 


(5.23) 
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6.0 ROTATION CONTROL PARAMETERS 

The attitude control matrix may be derived as from rotation 
equations as follows 


xi = ax^i 


X 2 = a 2 X 2 


Xi = Cl*! X 2 — c 2^2 


X! = [axCi] *1 X 2 = [a 2 c 2 ]X 2 t 6 * 1 ' 

The desired attitude condition at docking is that the two systems 
have approximately the same orientation in inertial space and, 
therefore. 


aiCi — 9-2^2 • 


( 6 . 2 ) 


If vehicle number one is the active control vehicle. 


aid — a 2 c 2 (c 1 )T 


(6.3) 


is the desired attitude matrix for vehicle number one. The 
control matrix becomes 


aiXi — ai d Xid 


Xi = [ (ai)Tai d ]X ld 


Xx = [ (ai) T a 2 c 2 (Ci)T]x ld 


C N = t (a x )T a 2 c 2 (C].)T] (Attitude Control Matrix) 


(6.4) 



That is, the C N matrix will become the unit matrix when the docking 
ports are aligned. By extracting the Hamilton Quaternion from the 
C N matrix a desired rotation can be found to align the docking 
ports as follows; 


r qi > 


r cos (a>/2) > 

<12 


sin (oo/2) cosa 

qa 


sin (<o/2) cosp 

<q4 , 


^sin (<o/2) cosy j 


where 


«o = 2 cos -1 (qi) , 


( 6 . 6 ) 


is the instantaneous angular control rotation required to align 
the docking ports . 


The eigen axis about which the 
must occur becomes 


< cosa > 

1 

/ \ 
q* 

cosp 

= 

qa 


sin (m/2) 


^cosy j 


< q<> 


controlled rotation of vehicle one 


(6.7) 


Hence, as a> approaches zero, the Cr matrix approaches the unit 
matrix and the desired attitude condition of Equation 6.4 can be 
reached. The vehicle one attitude GN&C functions must be designed 
to control the spacecraft to this desired attitude. 

The soft singularity indicated in Equation 6.7 will not be reached 
since a> will fall into the control attitude deadband before co «■ 0. 
A candidate commanded body axis rate vector may be estimated by 
selecting a scaler quantity n for the desired angular closing rate 
for the control quaternion and the desired body axis rotation rate 
becomes 
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n 


( 6 . 8 ) 


“id = 


cosa 



cosp 

a 

qa 

sin (co/2 ) 

^cosy j 


^ q-s j 


The candidate commanded body axis rotation rate change is 


A©i c = mi d “ «>if 


(6.9) 


where fi»i is the instantaneous body axis rotation rate vector for 
the active vehicle (i.e., vehicle one). 

The desired attitude rate constraint can also be derived by 
equating the two transformation matrices for docking port 
coordinates to the inertial coordinates as follows. 


(a 0 )T ai Ci = (a 0 )Ta 2 c 2 


( 6 . 10 ) 


and then demanding that the rates of the changed resulting 
matrices be identical, 


-W 0 (a Q ) T aiCi + (aJTaxWiC! = - W 0 (a 0 )fa 2 c 2 + (a Q )Ta 2 W 2 c 2 (6.11) 


and 


-a 0 W 0 (a Q ) Ta^i + a^Ci = -a 0 W 0 (a Q ) Ta 2 c 2 + a 2 W 2 c 2 . (6.12) 

Solving for the body axis rotation rate command matrix for vehicle 
one 


Wic = (a 0 )T{ (a 0 W 0 (a 0 )T) ( ai -a ld ) + a 2 W 2 c 2 (c x ) . 


(6.13) 
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The am matrix is the desired attitude matrix given by Equation 
6.3. Note that when the ai matrix approaches the a ld (desired 
attitude matrix) the steady state rotation rate command matrix for 
vehicle one becomes 


Wic = (ai) Ta 2 W 2 C2 (Ci) T . 


The candidate FCS could monitor the error signal, (ai~ai d ) from 
Equation 6.13 above and apply appropriate gains and limiters to 
control the attitude of vehicle one during docking. 
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7.0 CONCLUSIONS AND RECOMMENDATIONS 


Detailed analysis of the Automatic Rendezvous and Capture problem 
indicate three different regions of mathematical description for 
the GN&C algorithms; (1) multi-vehicle orbital mechanics to the 
rendezvous interface point, i.e., within 100 NM, (2) relative 
motion solutions (such as Clohessy-Wiltshire type) from the far 
field to the near field interface, i.e., within 1 NM and (3) close 
proximity motion - the near field motion where relative 
differences in the gravitational and orbit inertial effects can be 
neglected from the equations of motion. Limit boundaries to these 
regions can be precisely defined by further analysis and will also 
be a function of the tracking measurement accuracies and the 
computer resources available to the solution algorithms. 

The subjects of this paper analyze the regions (2) and (3) above 
and presents the derivation and discussion of the general case of 
gravitational perturbed relative motion. Mathematical deviations 
from the spherical gravity case of Clohessy-Wiltshire are 
presented in the analysis. Based upon the preliminary analysis of 
this work it is recommended that further efforts be used to assess 
the relative position and velocity differences in region (2) due 
to non-spherical effects (especially that of J 2 ) in the gravity 
perturbations. Future GN&C systems using relative motion 
measurements from either an onboard laser tracking systems or from 
the GPS will possibly be able to detect and predict these 
perturbations in the relative motion. 

The docking port relative position and velocity target vectors as 
outlined in this work couples the effects of the translation and 
rotation control activity. Based on these preliminary control 
parameters it is recommended that guidance and control systems 
functions couple translation and attitude control in region (3) 
for safe docking maneuvers. In region (2) the relative range 
between the docking port targets is assumed large enough so that 
translation and rotation G&C functions can be independent of one 
another . 
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SUGGESTED INTEGRABLE PARAMETER LIST 
APPENDIX A 


The following list contains the fundamental parameters necessary 
to solve for the motion and attitudes of two vehicles relative to 
the geodesic LVLH frame and the body axis frame, respectively. A 
numerical process to perform the first order solution of this set 
of simultaneous differential equations is required. This process 
will perform the integral 


ft+At 


x ( t +At ) = X ( t ) + 


x(t) dt . 


J t 


The process of Runge-Kutta-Gill is recommended, however, in some 
applications requiring high frequencies (i.e., small At) the Cauchy- 
Euler solution method is adequate for an accurate solution. 
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Suggested list 


x (t) 


x(t) 



* > 

1. 

2. 

r x oi 


3. 

4 . 1 

5. | 

1 

^ol 

Vehicle 

6. 

1 

Number " 

7- 1 

1 

One 

8. 

qi 


Vehicle 

Number 

Two 


9. 

10 . 
11 . 
12 . 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20. 
21 . 
22 . 

23. 

24. 

25. 

26. 


r ©! 


r *o 2 


r ?o 2 


r q 2 


©2 


(Position) 

LVLH Ref. 

(Velocity) 

LVLH Ref. 

(Quaternion) 
(Body axis to 
inertial) 

(Body Axis 
Rotation Rates) 


I Vol 

9o 

“Biqg o 0 V REL + Ti + An 
Wi 


r qi 


* Ji(Nl-WiIlfi»l) 

' -B 2 qg 0 UvREL + A i2 
k 42 


< 


J 


J 2 (N 2 -W 2 I 2 <b2) 


In the event that the absolute inertial state is required, the 
following additions can be made to the integrable parameter list 


x(t) 


*(t) 



27. 

1 


28. 

x ° r 

Inertial ’ 

29. . 

. 

State 

30. 1 



31. 

Vo 


32. J 

1 J 


where G 0 is the acceleration due to gravity at the point X Q , the 
geodesic coordinate center (see Figure 2-1) . 
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GRAVITATIONAL ACCELERATIONS 


APPENDIX B 


Using the Newtonian formulation of gravity the inertial 
gravitational acceleration vector in a spherical gravity field is 


G(x) = — jiR" 3 x 


(B.l) 


where x is the column vector for position from the gravitating 
mass center, that is. 


f 


x 


x 





J 


(B.2) 


To transform B.l into the LVLH coordinate frame use Equation 2.1 
in B.l, 


G(x) = - nR~ 3 a 0 X - fiR -3 x G 


(B . 3) 


and transform G(X) to the LVLH frame, 


(a Q ) T G (X) = G(X) = - tiR- 3 (X) - nR-3(aT Xo ) . 


(B.4) 


Since the LVLH coordinate axes are aligned with the Z-axis in the 
-R direction. Equation B.4 reduces to 


G (X) = - jiR 


-3 


^X ^ 


? 

V Z-R 0 j 


(B . 5) 
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The gravitational acceleration vector in a field perturbed by non- 
spherical gravity harmonics, we write. 


G (x) = - nR-3 (1+9?) x 


(B . 6) 


where I is the unit matrix an 9? is a diagonal matrix, 


9? = 




0 0 

9?y 0 

0 9? 


(B . 7 ) 


z J 


For example, using the J 2 harmonic only as shown in Reference 2, 
Pages 50 and 51, 
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(B . 8) 


Here, R E is the reference planetary equatorial radius and J 2 is the 
first zonal harmonic. For instance, the J 2 harmonic for the earth 
is -1 . 0826271xl0 -2 (unitless) and has the most pronounced effect on 
the motion of the low earth orbit spacecraft . 

Transforming G(x) in Equation B.6 into G(£) we have 


G (X) = -jiR-3 (1+9?) 


u 


(B. 9) 


where I is the unit matrix and 9? is given by the similarity 
transformation. 
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9t = a T SR a, 


(B. 10) 


and becomes a full matrix in the LVLH coordinate frame. 
J 2 harmonic only from Equation B.8 and reducing to LVLH 
Equation B.9 becomes. 


G (x) 


-jiR-3 


1+3/2 J 2 
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2 ' 


(x ) 

a T [K] a-5 
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1 R J 

> 

; 

^ 2^R 0 j 


where the matrix [K] is a constant matrix given by 


fl.O 0 0 


[K] 


0 1.0 0 . 


^ o o 3.0; 


Using the 
coordinates 


(B.ll) 


(B . 12) 
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PARTIAL DERIVATIVES 07 THE 
GRAVITATIONAL ACCELERATION VECTOR 

APPENDIX C 


The P x Matrix of Equation 2.15 for a spherical gravitational field 
is found by partial differentiation of Equation B.5 with respect 
to the x coordinates and is 
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(C.l) 


The Clohessy-Wiltshire relative motion equations result when X, y 
and 2 are very small relative to the radius vector magnitude R 
(thus R approaches R 0 ) and Equation C.l is approximated by 


( 1.0 0.0 0.0 


Px = -nR - 3 


0.0 


1.0 0.0 . 


(C . 2 ) 


^ 0.0 0.0 - 2.0 ) 


When the gravity field is perturbed due to non-spherical harmonics 
as in a real world situation the P x matrix is found by partial 
differentiation of Equation B.9 with respect to the X coordinates 
and can be written out in component form as follows : 


32 


3G* f _ 39? x _ 

= nR-3 l+* n + 3 — (l+*n) X 

3x l l 3x R 2 ) 


' 39?i2 X | 39?i3 X 

+ 3 — y + 3 — ( z-R 0 ) 

L 3X R 2 J l 3X R 2 ) 


3G X {_ 

r 3*n 

= jlR-3 *12 + 


3y V 

<3y 


-3 — (l+*n) X 

R 2 J 


[ a»!2 y 
+ — - 3 — *12 y 
lay R 2 J 


[3*i 3 y _ 

+ - 3 — *13 (2- R 0 ) 

lay R 2 J 


3G X f _ [ 3*11 (2-Ro) 

= ^R -3 *13 + 3 (H"*!!) X 

az 1 3z R 2 ; 


f a*i2 (2-R 0 ) 

+ 3 *12 y 

lay R 2 J 


(a*!3- (2-Ro) 

+ - 3 *13 (2 — R 0 ) 

. az R 2 


33 



X 


dG, 


f M 21 

2 


= jiR-3 

*21 + 

- 2 

9t 2 i 

9x 


^ dx 

R2 



^22 

X 

_ 1 


+ 

_ 

3 — ( l+95 22 ) 

y 


^ax 

R2 

J 



^*23 

2 > 



-l- 

3 

— *23 

(2-Ro) 


^ax 

R2 J 


J 


dGy 

f 

r d9t 2 i 

? _ ^ 

= >iR-3 

1 + SR22 + 

3 ■ 

*21 

ay 

< 

[dy 

R2 J 


f a* 22 y 

> 


+ 

— 3 — 

“ (I+9I22) 

y 


^ay R2 


r M 2 3 

^az 


3 




(C . 7) 


dGy 

az 



f 

'^*21 (2-Ro) > 

J»R-3 

*23 + 

- 3 <l+*2l) 



^ az R 2 > 


a 9 t 22 

^az 


(2-Ro) 


R2 


( 1+9I22) 


) 


^*23 

^az 


( 2— R 0 ) 

3 91 

R2 


23 


(2-Ro) 


(C . 8) 


34 


3G 2 

f 

3*31 

X 

= jiR-3 

*31 + 

- 

3 — *31 

3x 


^3x 

R 2 J 



( 3 SR 32 

X > 




+ 

— - 

3 — *32 

9 




^ ax 

R 2 J 





3*33 

X 

) 


) 

+ 

3 — (l+* 33 ) 

(Z-Ho) 


^3x 

R2 

) 


) 

3G 2 


f 3*31 

y 

\ 


- hR-3 

*32 + 

3 

— 

*31 

X 

ay 


3y 

R2 

J 



'3*32 

? _ ' 




+ 

- 

3 — *32 

y 




<3? 

R2 j 





' 3*33 

y 

> 

1 


+ 

- ; 

3 — (l+*33) 

(2-Ro) 


t,3y 

R 2 

j 

1 



(C . 9) 


(C.10) 



r 

fd* 31 (2-R 0 ) _ 

= hR'3 

1 + *33 + 

3 - 


*31 



^32 

R2 ; 


( 39132 (Z-Ro) _ 'l 



+ 

- 3 — 

*32 

9 



^32 

R2 J 




( 3*33 ( 2 - 

-Ro) 

) 

\ 

+ 

3 — 

(l+*33> 

(2-Ro) 


^32 

R2 

J 

/ 


(C.ll) 


35 



As X, y and Z approach zero (R approaches R 0 ) and Equations C.3 
through C.ll form the P x matrix as in Equation C.2 but for non- 
spherical case. 
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Typical values for the P x matrix in a 4x4 gravity field at 190 NM 
altitude (orbit inclination 28.5*) as a function of orbital 
longitude (L) are; 


For L = 0 . 0 * , 

P x = -HR- 3 

For L = 91. 914 # 
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(C.15) 
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